Steinmetz-Larter model
This example is found in the MatCont ecosystem.
The Steinmetz-Larter model is studied in the MatCont ecosystem because it is a simple example where a Chenciner bifurcation occurs.
\[\tag{E}\left\{\begin{array}{l} \dot{A}=-k_1 A B X-k_3 A B Y+k_7-k_{-7} A, \\ \dot{B}=-k_1 A B X-k_3 A B Y+k_8, \\ \dot{X}=k_1 A B X-2 k_2 X^2+2 k_3 A B Y-k_4 X+k_6, \\ \dot{Y}=-k_3 A B Y+2 k_2 X^2-k_5 Y, \end{array}\right.\]
This tutorial is also useful in that we show how to start periodic orbits continuation from solutions obtained from solving ODEs. Being a relatively advanced tutorial, we will not give too many details.
We start by coding the bifurcation problem.
using Revise
using Test, ForwardDiff, Parameters, Plots, LinearAlgebra
using BifurcationKit, Test
const BK = BifurcationKit
norminf(x) = norm(x, Inf)
function SL!(du, u, p, t = 0)
@unpack k1, k2, k3, k4, k5, k6, k7, k₋₇, k8 = p
A,B,X,Y = u
du[1]= -k1*A*B*X - k3*A*B*Y + k7 - k₋₇*A
du[2] = -k1*A*B*X - k3*A*B*Y + k8
du[3] = k1*A*B*X - 2*k2*X^2 + 2*k3*A*B*Y - k4*X+k6
du[4] = -k3*A*B*Y + 2*k2*X^2 - k5*Y
du
end
SL(u,p,t=0) = SL!(similar(u),u,p,t)
z0 = rand(4)
par_sl = (k1=0.1631021, k2=1250., k3=0.046875, k4=20., k5=1.104, k6=0.001, k₋₇=0.1175, k7=1.5, k8=0.75)
prob = BK.BifurcationProblem(SL, z0, par_sl, (@lens _.k8);)
argspo = (recordFromSolution = (x, p) -> begin
xtt = BK.getPeriodicOrbit(p.prob, x, set(getParams(p.prob), BK.getLens(p.prob), p.p))
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
period = getPeriod(p.prob, x, set(getParams(p.prob), BK.getLens(p.prob), p.p)))
end,
plotSolution = (x, p; k...) -> begin
xtt = BK.getPeriodicOrbit(p.prob, x, set(getParams(p.prob), BK.getLens(p.prob), p.p))
plot!(xtt.t, xtt[:,:]'; label = "", k...)
# plot!(xtt.t, xtt[2,:]; label = "y", k...)
# plot!(br; subplot = 1, putspecialptlegend = false)
end)We obtain some trajectories as seeds for computing periodic orbits.
using DifferentialEquations
alg_ode = Rodas5P()
prob_de = ODEProblem(SL!, z0, (0,136.), par_sl)
sol = solve(prob_de, alg_ode)
prob_de = ODEProblem(SL!, sol.u[end], (0,30.), sol.prob.p, reltol = 1e-11, abstol = 1e-13)
sol = solve(prob_de, alg_ode)
plot(sol)
We generate a shooting problem form the computed trajectories and continue the periodic orbits as function of $k_8$
probsh, cish = generateCIProblem( ShootingProblem(M=4), prob, prob_de, sol, 16.; reltol = 1e-10, abstol = 1e-12, parallel = true)
solpo = newton(probsh, cish, NewtonPar(verbose = true))
_sol = BK.getPeriodicOrbit(probsh, solpo.u, sol.prob.p)
plot(_sol.t, _sol[:,:]')
opts_br = ContinuationPar(pMin = 0., pMax = 20.0, ds = 0.002, dsmax = 0.05, nInversion = 8, detectBifurcation = 3, maxBisectionSteps = 25, nev = 4)
opts_po_cont = setproperties(opts_br, maxSteps = 60, saveEigenvectors = true, tolStability = 1e-3)
@set! opts_po_cont.newtonOptions.verbose = false
@set! opts_po_cont.newtonOptions.maxIter = 10
br_sh = continuation(deepcopy(probsh), cish, PALC(tangent = Bordered()), opts_po_cont;
verbosity = 3, plot = true,
callbackN = BK.cbMaxNorm(10),
argspo...) ┌─ Curve type: PeriodicOrbitCont
├─ Number of points: 61
├─ Type of vectors: Vector{Float64}
├─ Parameter k8 starts at 0.75, ends at 0.8117057297230886
├─ Algo: PALC
└─ Special points:
If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`
- # 1, ns at k8 ≈ +0.82521358 ∈ (+0.82521356, +0.82521358), |δp|=2e-08, [converged], δ = ( 2, 2), step = 21, eigenelements in eig[ 22], ind_ev = 2
- # 2, bp at k8 ≈ +0.84391360 ∈ (+0.84391360, +0.84391360), |δp|=5e-09, [converged], δ = (-1, 0), step = 32, eigenelements in eig[ 33], ind_ev = 2
- # 3, endpoint at k8 ≈ +0.81340853, step = 61
Curve of Fold points of periodic orbits
opts_posh_fold = ContinuationPar(br_sh.contparams, detectBifurcation = 2, maxSteps = 35, pMax = 1.9, plotEveryStep = 10, dsmax = 4e-2, ds = 1e-2)
@set! opts_posh_fold.newtonOptions.tol = 1e-12
@set! opts_posh_fold.newtonOptions.verbose = true
fold_po_sh = @time continuation(br_sh, 2, (@lens _.k7), opts_posh_fold;
verbosity = 3, plot = true,
detectCodim2Bifurcation = 2,
startWithEigen = false,
usehessian = false,
jacobian_ma = :minaug,
normN = norminf,
callbackN = BK.cbMaxNorm(1e1),
bothside = false,
bdlinsolver = BorderingBLS(solver = DefaultLS(), checkPrecision = false),
)
plot(fold_po_sh)
Curve of NS points of periodic orbits
opts_posh_ns = ContinuationPar(br_sh.contparams, detectBifurcation = 1, maxSteps = 35, pMax = 1.9, plotEveryStep = 10, dsmax = 4e-2, ds = 1e-2)
@set! opts_posh_ns.newtonOptions.tol = 1e-12
ns_po_sh = continuation(br_sh, 1, (@lens _.k7), opts_posh_ns;
verbosity = 3, plot = true,
detectCodim2Bifurcation = 2,
startWithEigen = false,
usehessian = false,
jacobian_ma = :minaug,
normN = norminf,
callbackN = BK.cbMaxNorm(1e1),
bothside = false,
) ┌─ Curve type: NSPeriodicOrbitCont
├─ Number of points: 36
├─ Type of vectors: BorderedArray{Vector{Float64}, Vector{Float64}}
├─ Parameter k7 starts at 1.5, ends at 1.8383023734828383
├─ Algo: PALC
└─ Special points:
If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`
- # 1, ch at k7 ≈ +1.74714606 ∈ (+1.74714588, +1.74714606), |δp|=2e-07, [converged], δ = ( 0, 0), step = 27, eigenelements in eig[ 28], ind_ev = 0
- # 2, endpoint at k7 ≈ +1.84894792, step = 36
plot(ns_po_sh, fold_po_sh, branchlabel = ["NS","SN"])